function [f,f_n,Aineq_d,Aineq_i,bineq,lb,ub,ctype, Ind] = generate_cplex_inputs(bsperm,g,X,bsi,k,gn,Xn)
% This function generates inputs required for the CPLEX solver

% Input: 
% (1) bsperm: bootstrap permutation index 
% (2) g: treatment effect calculated using IPW 
% (3) X: set of covariates 
% (4) bsi: bootstrap index 
% (5) k: number of covariates 
% (6) gn: original sequence of g 
% (7) Xn: original sequence of Xn

% Output 
% (1) f: objective function coefficients
% (2) f_n: reverse objective function coefficients
% (3) Aineq_d: constraint matrix in decreasing order 
% (4) Aineq_i: constraint matrix in increasing order
% (5) bineq: constraint right hand side
% (6) lb: lower bound on solutions (beta and z)
% (7) ub: upper bound on solutions (beta and z)
% (8) ctype: variable type 
% (9) Ind: sort index of households by X in original sample or bootstrap
% sample

%% Set up cplex inputs 
    gr = g(bsperm,:);
    Xr = X(bsperm,:);
    % Subtract Wn from Wbs
    if (bsi>0)
        gr = [gr; -gn]; 
        Xr = [Xr; Xn];
    end
    % Compress data with identical X's
    [Xr, Ind] = sortrows(Xr);
    gr = gr(Ind);
    Xu = unique(Xr,'rows');
    nu  = size(Xu,1); % number of unique covariate vectors
    gu = zeros(nu,1);
    jj = 1;
    for j = 1:size(Xr,1)
        if ~(sum(Xr(j,:)~=Xu(jj,:)))
            gu(jj) = gu(jj) + gr(j);
        else
            jj = jj+1;
            gu(jj) = gu(jj) + gr(j);
        end
    end

    % Add explicit monotonicity constraints
    % Xu is ordered by increasing covariate (Xu(:,2))
    % and then by increasing usage (Xu(:,5))
    % For each level of covariate we impose treatment set inclusion
    sameinc = (Xu(1:end-1,2)==Xu(2:end,2));
    % decreasing in pre-treatment consumption
    Mineq_d = [diag(sameinc) zeros(nu-1,1)] + [zeros(nu-1,1) diag(-sameinc)];
    % increasing in pre-treatment consumption
    Mineq_i = [diag(-sameinc) zeros(nu-1,1)] + [zeros(nu-1,1) diag(sameinc)];

    f = [zeros(k,1); gu]; % objective function coefficients
    f_n = [zeros(k,1); -gu]; % Reverse objective function coefficients
    B = 1; % bounds on coefficients
    C = B*sum(abs(Xu),2); % maximum values of x'beta
    minmargin = max(1,C)*(1e-8); % prevent non-integer numbers the integrality constraint of integers from being counted as integers
    Aineq_d = [[-Xu diag(C)]; [Xu -diag(C)]; [zeros(nu-1,k) Mineq_d]]; 
    Aineq_i = [[-Xu diag(C)]; [Xu -diag(C)]; [zeros(nu-1,k) Mineq_i]];
    bineq = [[C-minmargin];[-minmargin]; minmargin(1:nu-1,:)];
    lb = [-B*ones(k,1); zeros(nu,1)];
    ub = [ B*ones(k,1);  ones(nu,1)];
    
    % Variable type string
    ctype = strcat(repmat('C',1,k),repmat('B',1,nu));

end 